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Abstract 

Given a marked renewal point process (assuming that the marks are i.i.d.) we say that an unbounded 
region is stable if it contains finitely many points of the point process with probability one. In this 
paper we provide algorithms that allow to sample these finitely many points efficiently We explain 
how exact simulation of the steady-state measure valued state descriptor of the infinite server 
queue follows as a simple corollary of our algorithms. We provide numerical evidence supporting 
that our algorithms are not only theoretically sound but also practical. Finally having simulation 
optimization in mind, we also apply our results to gradient estimation of steady-state performance 
measures. 

1 Introduction 

Let N = {N {t) : t G (— oc,oc)} be a two sided time stationary renewal point process. We write 
{An : n G Zq} for the times at which the process N jumps, where Zq = ^\{0} denotes the set of 
integers removing zero, and with Ai > > A-\. For simplicity we assume that A^ < ^n+i ioi 
every n. Further, we define = ^n+i — ^n- 

Now let {Vn : n G Zq} be a sequence of independent and identically distributed (i.i.d.) random 
variables (r.v.'s) which are independent of the process N. Define = (An^Vn) and consider 
the marked point process = {Z^ : n G Zq} which forms a subset of M?. We say that a (Borel 
measurable) set B is stable if \M nB\ < oo almost surely (where \C\ is used to denote the cardinality 
of the set C). 

Under natural assumptions on the inter-arrival times underlying N and on the distribution of 
the VnS (stated in Section [2]) we propose and study a class of algorithms that allow to sample 
exactly (i.e. without any bias) a realization of the set MdB for a large class of unbounded, stable 
sets B. 

Our approach builds on algorithms that are fully developed and studied in [4J. As an application 
of the class of algorithms that we study here, we provide a procedure that allows to sample from 
the steady-state measure valued descriptor of an infinite server queue without any bias (i.e. exact 
simulation). Such a procedure, for instance, is obtained by considering the particular case in which 
B takes the form B = {{t^v) : v > \t\,t <0}. Given that point processes constitute a natural way 
of constructing queueing models in great generality, we believe that the class of algorithms that we 
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propose here have the potential to be apphcable to the design of exact samphng algorithms of more 
general queueing models. This is a research avenue that we plan to investigate in the future. 

We argue empirically that it is cheaper to run our exact sampling procedure to fully delete the 
initial bias than it is to do a burn-in period that reduces the bias to a reasonable size, say 5%, when 
talking about, for instance, the steady-state queue length. 

Finally, we apply our exact sampling algorithms for infinite server queues to perform steady-state 
sensitivity analysis. For instance, we consider quantities such as the derivative of the steady-state 
average remaining service time with respect to the arrival rate or service rate. These quantities are 
of great interests in stochastic optimization via simulation. 

So, in summary, our contributions are as follows: 

i) We provide the first exact sampling algorithm for stationary marked renewal processes on 
unbounded and stable sets, see Section [2j 

ii) As a corollary of i) we explain how to obtain an exact sampling algorithm for the steady-state 
measure valued descriptor of the infinite server queue. We also show empirically that this algorithm 
is practical in the sense of being both easy to code and fast to run, see Section [3) 

iii) Finally, we provide new procedures for the sensitivity analysis of steady-state performance 
measures of the infinite server queue, see Section [4j 

Relevant literature 

Following the seminal work by [TT], several exact sampling algorithms have been developed, partic- 
ularly for spatial point processes. [9J and [lOj developed algorithms and analytical tools based on 
so-called Dominated Coupling From the Past (DCFP). DCFP is based on the idea of introducing a 
stationary dominating process that is simulatable. Compared to our method, firstly they use spatial 
birth and death processes (generally of poisson type) as the coupled dominating processes. This 
would limit the target distribution to be absolutely continuous with respect to the Poisson measure. 
Secondly the number of steps simulated in the naive DCFP grows exponentially with the system 
scale (i.e. arrival rate in the infinite server queue setting); see Proposition 1 in [3J for a detailed 
proof. Although several modifications have been proposed, still the number of steps involved in 
these backward construction appears to be significantly large, especially when sampling in infinite 
volume regions [7]; see Section 7 in [3] for empirical comparisons. 

Our method is based on a construction that is being used in [5j and [4J ; see also [6J for related 
ideas. The method involves the technique of simulating the maximum of a negative drift random 
walk and the last passage time of independent and identically distributed random variables to an 
increasing boundary. As shown in [4J the complexity of our algorithm scales graciously as the system 
scale grows. 

2 Sampling form stable unbounded regions 

We start by discussing the assumptions behind our development. 
Assumptions: 

Al) Assume that E\Vn\^^^ < oo for some > 0, we also write F {•) = P {Vn < •) for the 
cumulative distribution function (CDF) of Vn and put F (•) = 1 — F (•) for the tail CDF. 

A2) We assume that F (•) is known and easily accessible either in closed form or via efficient 
numerical procedures. Moreover, we can simulate Vn conditional on E [a, b] with P (Vn E [a, b]) > 
0. Finally we can find u{k) such that u{k) > PdVi]^/*^ > v)dv and u{k) ^ as /c ^ oc. 
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A3) Recall that = ^n+i — > 0. Define (9) = logE'exp {OX^) and assume that there 
exists 6 > such that t(j (6) < oo. Finally, let us write fi = EX^. 

A4) Define G {•) = P {Xn < •) and G (•) = 1 - G (•). Suppose that G (•) is known and that it is 
possible to simulate from Geq (•) := G (t) dt. Moreover, let Gq (•) = E exp {0Xn — tjj (6)) I {Xn < •) 

be the associated exponentially tilted distribution with parameter 9 for (9) < oo. We assume that 
we can simulate from Go {•). 

Consider the class of sets B cM? that are Borel measurable and such that 

BcCa^{{t,v) : \v\ > \tn. 

Our goal in this section is to develop an algorithm that allows to sample without bias the random 
set M nCa^ and therefore M D B. We will discuss extensions that follow immediately from our 
formulation at the end of this section. Figure [T] illustrates the different shapes that the set Ca can 
take depending on the values of a > 0. 



a=\ a>l 0<a<l 

Figure 1: The area of Ca- The horizontal axis corresponds to the t coordinate while the vertical 
axis represents the v coordinate 

We now proceed to explain our construction. As the stationary renewal point process is time 
reversible, starting at the distribution of the forward process {Z^ : n > 0} and the backward 
process {Z^ : n < 0} are the same. In what follows we limit our discussion to the construction of 
the forward process and the simulation of the backward process is completely analogous. 

We follow the idea in Let e G (0,/i). Consider any random time finite with probability 
one but large enough such that 

An+i > n{fi - e) and |K+i| < {n{fi - e))"^ 

for all n> n. 

If such random time is well defined, we only need to simulate the stationary process up to 
to get a sample from the unbounded region. 

Proposition 1. The random time n defined above exists and it is finite with probability one. 

Proof. By Chebyshev's inequality, 

P{An+i < n{fi - e)) < E[eMO{n{fi - e) - An+i))) < exp(-n(-0(^ - e) - ^{-9))) 

for any > 0. 
Let 

/(-6) = max{-6>(/i - e) - ^(-6>)} 



As V^(0) = 0, V^'(O) = fi and f'{0) = Var{X) > 0, /(-e) > 0. Then 

P{An+i < n{fi - e)) < exp(-n/(-6)) 

and 

t P(^*^ <n(,-e))< , !Tx;(-7(-')) ~ 



n=l 



By Borel-Cantelli lemma, {A^+i > — e)} eventually almost surely. 
Similarly and independently we have 

^P(|K+i|>(n(/x-6)n = ^P(|yi|^/^>n(/x-6))<- / Pm\'/^>u)du 

n=l n=l ^ " ^ "^0 



< OO 



Thus, again by Borel-Cantelli lemma, ^ (^(/^ — ^))^} eventually almost surely. Therefore, 

P{k < oo) = 1 □ 

As {An : n > 1} and {1^ : n > 1} are independent of each other, we consider the following 
construction. Let i^{A) be a random time satisfying that A^+i > — e) for n > hc{A), and 
be a random time satisfying that Vn+i < n{fi — e) for n > i^{V). Clearly Hi{A) and i^(V) are not 
stopping times and this makes the simulation of these times challenging. However, we will explain 
how to sample these times and then we can set k = max{K:(A), t^{V)}. Our construction will allow 
us to simulate {A^ : n > 1} and {Vn : n> 1} separately. 



2.1 Simulation of {A^ - 1 < k < max{n, k (A)} + 1} 

In this subsection we will introduce a method to simulate i^{A) together with {A^ : A: > 1}. 

First, define Ai according to the distribution Geq (•)• Sampling Ai can be done according to 
A4). 

Now, observe that A^+i = Ai + Xi + ... + and define 

n 

Sn = n(/i - e) - {An+i - Ai) = ^ Yi, 

1=1 

where Yi = [ji — e) — Xi. Note that the l^'s are i.i.d. with EYi = — e. If we set = 0, then 
{Sn : n > 0} is a random walk with negative drift. We are interested in sampling up to the last 
time n at which > 0. 

We define the following sequence of random times: 

Ai = 0, Ti = inf{n > Ai : - > 0}, 

and for j > 2 

A^- = inf{n > r^_il{r^_i < oo} V A^_i : Sn < 0}, 
r^- = inf{n > A^- : Sn - S^. > 0}. 

Now, let 7 = infjj > 1 -.Tj = oo} and note that A^+i = A^ and that < for n > A^, which 
in particular implies that A^+i > n(/i — e) for n > Aj. Therefore, we have that = k (A). 
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In what follows we will explain how to simulate the Aj's and Fj's sequentially and jointly with 
the underlying random walk until time A^. One important observation is that for every j > 1, 
Aj < oc almost surely by the strong law of large numbers. 

Let us write — ^2, Yn} ioi the cr-field generated by the Y^'s up to time n. Let ^ > 

and define 

:= inf{n > : 5^ > C}, 
then by the strong Markov property we have that for j < 7, 

P(r, = oo\Ta,) = P(r, = oo\Sa,) = P{n = oc) > 0, 

where we use P {•) to denote the nominal probability measure under which Sq = 0. 
It is important then to note that 



P = k) = P {To < 00)''-^ P {To 



OC 



for A: > 1. In other words, 7 is geometrically distributed. The procedure that we have in mind is 
to simulate in time intervals, and the number of time intervals is precisely 7. 

Let iPy{^) — log E ex.p{9Yi). As the moment generating function of Xi is finite in a neighborhood 
of the zero, '0y (•) is also finite in a neighborhood of zero and EYi = '0y(O) = — e, Var(y^) = '^y(O) > 
0. Then by the convexity of ipri'), one can always select e > sufficiently small so that there exists 
77 > with tl^riv) — ipyi^) ^ ^- "^^^ ^^^^ ^ allows us to define a new measure based on 
exponential tilting so that 

rIP 

^(y.) = exp(r?r.)- 

Moreover, under P^, Sn is random walk with positive drift equal to '0y (77) ([Ij P. 365). Therefore 
P^(To < oc) = 1 and P(To < oc) = Erj{ex.p{-r]STo))- More generally, Pry(T^ < oc) = 1 and 

q{0 Pm < oc) = E^{exp{-r]ST,)) 

for each ^ > 0. Based on the above analysis we now introduce a convenient representation to 
simulate a Bernoulli random variable J (^) with parameter q (^) namely, 

J{0 = I{U<exp{-r^ST,)). (1) 

where U is & uniform random variable independent of everything else under P^. 

Identity Q provides the basis for an implementable algorithm to simulate a Bernoulli with 
success probability q(^). Sampling {Si, Stq} conditional on To < oo, as we shall explain now, 
corresponds to basically the same procedure. First, let us write P*{-) — P{-\To < oo). The following 
result provides an expression for the likelihood ratio between P* and P^. 

Lemma 1. We have that 

^^*-{Si St ) = '^^P(~^'^To) < 1 



Proof. 

P{Si(EHi,...,Sto(^Hto\To<oo 



dP^' p(ro<oo) -p(ro<oo)' 

P{Si(eHu...,Sto^Ht,,To<^) 



P{To < oo) 

Er,[exp{-riST,)I{So £ Ho,...,Sto £ Htq)] 
P{To < ^) 

□ 
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The previous lemma provides the basis for a simple acceptance / rejection procedure to simulate 
{Si, Stq} conditional on Tq < oc. More precisely, we propose (Si, Stq) from (•). Then one 
generates a uniform random variable U independent of everything else and accept the proposal if 

1 dP* ~ ~ 

^ ^ l/P{n < oo) " ^(5i,-,5to) = exp(-,5To). 

This criterion coincides with J (0) according to ([T]). So, the procedure above simultaneously obtains 
both a Bernoulli r.v. J (0) with parameter g'(O), and the corresponding path {S^i, S^Tq} condi- 
tional on To < oc. 

Algorithm 1 (Outputs (S^o, . , 5a^)) 
Step 0. Set K = 0, and Sq = 

Step 1. Simulate (Si, ...,Sto) from P^ and compute J := J (0) according to 

Step 2. If J = 1, then let Sx+j = Sj for j = 1, ...,To and update K ^ — K + Tq. Then, go back to 
Step 1. 

Otherwise, J = (i.e. = K), stop and output (S^o, Sk) 

Remark: It has been proved in [3] that the expected number of times we need to repeat Step 1 
does not change with the system scale (i.e. the arrival rate). 

We noted earlier that = (A) and Algorithm 1 together with the initial procedure to 
sample Ai allows us to simulate (^j+i : < j < (A)), and we know that A^+i > n(fi — e) for 
n > Hi (A). We need to simulate A^+i for n < = max{/^ (A) ,hc(V)}, and i^(V) is independent 
of hi (A). So, there might be cases for which we will have to sample A^+i for n > hi (A). Since 
An^i = Ai — Sn + n{fi — e) it suffices to explain how to simulate Sn for n > A^. In turn, it suffices to 
explain how to simulate {Sn n > 0) with S^o = conditional on Tq = oc. We will once again apply 
an acceptance / rejection procedure but this time we will use the original (nominal) distribution as 
the proposal distribution. Define 

p'(-) = P(-|ro = oo). 

The following result provides an expression for the likelihood ratio between P' and P. 
Lemma 2. We have that 



Proof. 

P{SieHi,....,SieHi\To = ^ 



dP' ' ' P(ro = oo) -p(ro = oo) 

P{SieHi,...SieHi,To = (x>) 



Pin = oo) 

E[IiSi G Hi,..., Si G Hi)IiTo > l)P{To = (x>\So, Si)] 



P{To = oo) 

The result then follows from the strong Markov property and homogeneity of the random walk. □ 
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We are in good shape now to apply acceptance / rejection to sample from P^ The previous 
lemma indicates that to sample {Sq, 5^} given Tq = oo we can propose from the original (nomi- 
nal) distribution and accept with probability q{—Si) as long as Sj < for all < j < L So, in order 
to perform the acceptance test we need to sample a Bernoulli with parameter q{—Si), but this is 
easily done using identity ([T]). Thus we obtain the following procedure. 

Algorithm 2 (Given n > outputs {^i, ^2, -4^ax{n,^(A)}+i}) 
Step 1. Run Algorithm 1 and obtain {S^o, S'l, Sk}- 

Step 2. U K = n (A) > n, jump to Step 6. Otherwise, K < let I = n — K > 1. 

Step 3. Simulate {S^o, S^i, 5/} from the original (nominal) distribution with Sq = 0. 

Step 4. If Sj < for all < j < / then sample a Bernoulli J{—Si) with parameter q{—Si) using 
and continue to Step 5. Otherwise (i.e. Sj > for some 1 < j < /) go back to Step 3. 

Step 5. If J{—Si) = 1, go back to Step 3. Otherwise, J{—Si) — 0, let Sx+i = Sk + Si for z = 1, 2, / 

Step 6. Let m = max{n, /^(A)}. Simulate Ai with CDF Geq{-) = fi~^ G{t)dt. Set A^+i = Ai - 
Sn + ^(m ~ ^) for n = 1, m. Output {Ai, A^+i}. 

2.2 Simulation of {Vn : I < n < i^{V) + 1} 

In this section we will introduce a method to simulate i^{V) together with the {Vn : n > 1}. 

Let p{n) = P{\Vi\ > {n{ii - e))^). We define Tq = and = inf{n > Ti_i : |14+i| > (n(/x - 
e))^} for i = 1, 2, .... We also define two independent sequences of random variables, {V^+i : n > 1}, 
and {V^+i : n > 1} as follows. The elements in each sequence are i.i.d., is distributed as Vn+i 
conditional on |V^+i| > {n{ji — e))*^, and V^+i follows the distribution of T4i+i conditional on 
^ (^(m ~ ^))^- We simulate Vi following its nominal distribution independent of everything 

else. 

Let (J — inf{z > : = oo}. Then l^i+i ^ (^(m ~ ^))^ for n > T^r-i + 1. We next introduce a 
method to sample Ti, T2, ... sequentially and jointly with the V^s up until Tcr_i. 

The following lemma provides the basis to guarantee the termination of our procedure. 

Lemma 3. If E\Vi\^/'' < 00, then 

00 

P(Ti = oo) = - p{i)) > exp(-2£;|yi| V-/(^ - e)) > 0, 

i=l 

consequently Ea < exp(2£;|y | V«/(^ - e)) < oc. 

Remark: The bound on Ea can be improved. This improvement is important for the theoretical 
asymptotic analysis of GI/GI/oc application, see 

Proof. 

00 

p(Ti=oo)=n(i-^'H) ^ 

n=l 

> 



YleM-Mn)) 



n=l 



exp(-- 



f^-^ Jo 



poo 

Jo 



> v)dv) = exp(— 



2E\Vi\^/'' 
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For i = 2, 3, ... conditional on T(i — I) = k: 

P(Ti = oo|Ti_i = fc) = (1 -p(n)) > exp( > exp( ) 

n=k+l ^ ^ 

Thus a is stochasticahy dominated by a geometric random variable with parameter p = exp(— 2£'| Vi |^/^/ (/i— 
e)), the result then follows. 

□ 

Notice that 

n (1 > m^ = oo|T,_i = A:) > n (1 X exp(--^^ ^) (2) 

for / > A: + 1. 

Thus if we are simulating / ^ Bernoulli(r^) with := P(T^ = (X)|T^_i), then with probability 
one we can check whether U < P{Ti = oc|T^_i) for U ^ Unif[0, 1] by making / sufficiently large 
without calculating the infinite product in the definition of P{Ti = oc|Ti_i). 
On the other hand, if we define flj^ill ~ pU)) 1? ^^en 

P^T IT ^ ^ r I 1 

P(Ti = n|Ti < oo) = p(-^^<^) < ^'(")p(T,<oo)- 

Consider a random variable A'' with the following probability density function 

P{N = n) = cp{n) 

for n = 1,2,..., where c = (E^iP(^))"^- Then P(Ti = n|Ti < oo)/P(A^ = n) < l/(cP(Ti < 
oo)). 

So we can simulate Ti given Ti < oc using acceptance / rejection with N as the proposal 
random variable. Generalizing the idea to T^, we can obtain the following algorithm 

Algorithm 3 (Given T^_i = A:, outputs conditional on < oo) 

Step 1. Let c = {J2^=k^i P{^))~^ • Simulate N with probability density function P(N — n) — cp{n) 
forn = A: + 1, A: + 2, ... 

Step 2. Simulate U - Unif[0, 1] independently, li U < Y{^S^j^i{l - p{j)) , set = TV and stop. 
Otherwise go back to Step 1 

We conclude this section with our procedure to simulate {Vi^ V2, •'•V^(v)+i}- 

Algorithm 4 (Outputs {Vi, ^2, 
Step 0. Set To = 0, z = 1. Simulate Vi from its nominal distribution. 
Step 1. Simulate / ^ Bernoulli(r^) with := P(T^ = oc|T^_i) (see Q). 

Step 2. If / = 1, set i^(V) = T^_i + 1. Simulate V^(y^^i by sampling from V^(v)^i and stop. Otherwise 
/ = 0, sample conditional on < oc and the value of T^_i using Algorithm 3. Simulate 
the process between T^_i + 2 and + 1 by sampling from Vn for T^_i + 2 < n < and Vn 
for n = + 1. Set i = z + 1 and then go back to Step 1. 
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3 Application to the infinite server queue 

As a direct application of the ideas discussed in the previous section we study steady-state simulation 
for the infinite server queue. The following diagram indicates how to construct the steady-state 
measure valued descriptor assuming that we can sample all the points inside the set 

C = {{t,v) :v>\t\,t<0}. 

Let Q{t,y) denote the number of people in the system at time t with residual service time strictly 
greater than y and E(t) denote the time elapsed since the previous arrival at time t (i.e. E (•) is the 
age process associated with N {•)). Figure [2] below depicts the region C. Every point in \A4 nC| is 
projected to the vertical line at time zero by drawing a —45^ line. The final position in the vertical 
line if positive, represents the corresponding remaining service time. Since the underlying point 
process is time stationary, the whole configuration of points obtained by this procedure at time zero 
is a snap shot of the steady-state distribution of the infinite server queue. 



\ \ Q(t,y)' 










y 








Figure 2: The points lies in the shaded area correspond to people who are still in the system at 
time with remaining service time greater than y 

3.1 Algorithm for the Infinite Server Queue 

As depicted in Figure [2] after projecting into the vertical line at t = 0, we obtain the stationary 
remaining service requirements of the customers at time zero. We shall use i?2, i?g(o,o) 
denote the remaining service times. The labeling is arbitrary although we will assign smaller indexes 
to customers that have spent less time in the system. Our algorithm proceeds as follows. 

Algorithm 5 (Outputs i?2, i?Q(o,o)} and £"(0)) 
Step 1. Use Algorithm 4 to simulate the {Vn, 1 < ^ < ^(^) + !}• 
Step 2. Use Algorithm 2 to simulate the {Ai^A2^ •••5 ^max{At(y),At(A)}+i}- 

Step 3. Set hc = max(/^(y), /^(A)). If > /^(V^), simulate Vn by sampling from for n = i^{V) + 
2,...,/^+!. 

Step 4. Set g = 0, z = and repeat the following procedure until i — Kj'. 
set i — i^l;ifVi> Ai, set q — q^l and Rq — Vi — Ai. 
Output i?2, •••^g} and Ai. 
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3.2 Empirical Performance 

Let y = {y(t) : t > 0} be a continuous time Markov process on the state space Q and / is a 
real-valued function defined on fi. The ergodic theorem guarantees in great generality (assuming a 
unique stationary distribution tt (•)) that 

J f f{Y{s))ds^ [ f{y)7:{dy) 
^ Jo Jn 

as t ^ oo almost surely for every positive, measurable function / (•). In the setting of the infinite 
server queue such a stationary distribution exists if EVn < oc and EX^ < oc. The most natural 
estimator for E^f{Y) := j^f{y)Ti{dy) is therefore 

$(t,y(0)) ■.= \ jj{Y{s))ds, 

where 1^(0) is the initial state. The estimator $ (t, 1^(0)) is generally biased unless Y (0) is sampled 
from the stationary distribution tt (•) ([2J P. 97). Our algorithm has the obvious advantage of 
removing the initial transient. 

In what follows we conduct some simulation experiment to evaluate the practical performance 
of our algorithm. The idea is to fix a reasonable tolerance error, say 10%, for a given performance 
measure. Then we want to empirically find how large a burn-in period one would need in practice 
to reduce the initial transient bias to about 10%. In order to effectively quantify the error we select 
a class of systems for which tt (•) can be explicitly evaluated. 

We consider an infinite server queue with Poisson arrivals and Lognormal service times. As we 
are interested in the efficiency of our algorithm for relatively large systems, we set the arrival rate 
A = 100 and the service time Vn ^ Lognormal(— 0.25, 0.5) (i.e. Vn has the same distribution as 
exp (—.25 + .5 X A^(0, 1)), where (0, 1) denotes a standard Gaussian random variable). 

Let Y{t) — {Q{t, •), E(t)) G P[0, oc) x R+, then Y(t) is a Markovian measure valued descriptor 
of the infinite server queue (of course in the Poisson arrival case one does not need to keep track of 

We first compare the performance of our algorithm to the burn-in period defined as the period 
needed to reduce the initial transient as indicated earlier. Let f{Y{t)) = (5(t,0), i.e. the number 
of people in the system at time t. We measure the computation effort of the algorithm in terms 
of the number of arrivals (we call this the number of steps) simulated. Given e > we let n(e) 
denote the minimum number of steps required so that \E^{An(^^^^ (0, 0)) — £^7rQ(0, 0)\/E^Q{0, 0) < e, 
where (0, 0) denotes a system that starts empty with E{0) = (recall that E{-) is the age process 
associated with A^(-), i.e. when £"(0) = x, Ai is distributed as Xn conditional on Xn > x). Table [l] 
shows the relation between e and n(e), obtained empirically based on the average of 10^ independent 
replications 



Table 1: Bias of $(S'^(,)) 



€ 


n{e) 


computer time (s) 


10.26% 


6 X 10^ 


0.0310 


5.71% 


1 X 10^ 


0.0382 


1.17% 


5 X 10^ 


0.1367 
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Compared to the results in Table [TJ our algorithm is unbiased. The average number of steps 
involved is n = 592.6369 based on the average of 10^ independent replications and the average 
computer time needed for a single replication is 0.0249 s. 

In addition, in Table[2]we compare the performance of the estimators (0, 0)) and ^(A^/, ((5(0, •), Ai)), 

where (5(0, •) and Ai are sampled according to Algorithm 5. n and are calibrated so that the 
computation budget is basically the same in both estimators. Under our procedure, Ehc, the average 
number of arrivals required to terminate is approximately equal to 600. So for instance, the first 
row in Table 2 corresponds to n = 10"^. This means that ^ 9.4 x 10^ = 10"^ — 600. The true 
value of Etj;Q{0^ 0) is 88.2497. The sample mean and sample standard deviation are calculated using 
the method of Batch means. The result in Table [2] shows that our mixed method performs better 
than the batch means with relatively small computation budget, while with large budget, the two 
methods are about the same. 

Table 2: Simulation result with different initial states. 





(<^,0) 


(Q(o,-) 




n 


Sample Mean 


Sample Std 


Sample Mean 


Sample Std 


1 X 10"^ 


86.1274 


1.0104 


88.1713 


0.6018 


5 X 10"^ 


89.0893 


0.4587 


88.2956 


0.3770 


1 X lO'^ 


88.5151 


0.3531 


88.1270 


0.2976 


5 X lO'^ 


88.3022 


0.1481 


88.3581 


0.1402 



4 Application to sensitivity analysis of the infinite server queue 

In this section, we apply our algorithm to sensitivity analysis of the infinite server queue. We 
consider a sequence of systems indexed by (A, z/), A > 0, z/ > 0. Given (A, z/), the interarrival times 
are multiplied by 1/A, obtaining X^/A for all n, and the service times are multiplied by thus 
we have Vn/i^ for all n. We assume that EVn < oc and EXn < oc. We will use the notation 
Qx,iy (•) to denote the infinite server queue descriptor for the (A, z^)-system. Our strategy rests on 
the application of Infinitesimal Perturbation Analysis (IPA), see for instance [8J P. 386. We assume 
here that the interarrival times have a continuous distribution. 

We illustrate the methodology by computing the sensitivity of the steady-state average remain- 
ing service time, which we denote by Et^R{\^ z/); namely, 

E^R{\ v) = E^ \ yQx^, (0, dy) . 

We also consider 

E^R^{\ v) = E^{\^i{y > : Qx,. (0, y) = 0}), 

in words, the steady-state maximum remaining service time. In order to apply IPA we need to 
define a few quantities. 

First, let us define S(A, u) to be the average elapsed service time of the customers that are 
present at time zero (given the construction of the stationary process {Qx^^ (t, •) : t G (— oo, oo)}, 
see Figure 2). That is. 
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Likewise, define ^ (A, u) as the average of the total service requirement of the customers that are 
present at time zero, namely 



Next, we define (A, v) as the elapsed service time of the customer with the maximum remaining 
service time at time zero and ^(^^(A, z/) as his total service time requirement. Specifically, if we 
let m = argmaxjn : V^/z^ — |A^|/A} then 

S(-) (A, I.) = ^ and F(°°)(A, z/) = ^ 
A V 

We then obtain the following representation for the derivatives of Ej^R{X^i/) and Ej^R'^{X,i/) 
with respect to A and u. 

Lemma 4. We have that 

^E^R{\v) = jE^E{\,i^) and ^E^R{\v) = -^E^V{\v). 

ii) 

^E^R^{X,u) = \e^E(^\X,u) and ^E^R^{X,u) = -^E^V^^\X,u) 

OA A OV V 

Proof. We only give a proof of part i) here as the proof of part ii) is entirely analogous. 

Let Rn denote the remaining service time of the nth customer at time zero and Vn as his total 

service time requirement, then Rn < Vn- Thus if EVn < cxo, we have 

E^R{X, v) < oo 

for any A > 0,z/ > 0. 

For a fixed sample path uo constructed backward in time, let i2ri(A, z^, cj), n < 0, denote the re- 
maining service time of customer n (counting backward in time) at time in system (A, v). Then 
Rn{\T^,uj) = {Vn{uj)/iy-\An{uj)\/Xy and 

i^n(A + h, V, UJ) - Rn{X, V, CJ) ^ |^n(^) | ^ p K(^) ^ |^n(^)| ^ 

/i^o h V ~ X 



i^n(A, iy + h,Uj) - RnjX, V, Uj) ^ Vn{^) ^^ Vn{^) ^ |^nM| 

/i-^o h v'^ V ~ X 



Thus the derivative -^R{X, v) and -^R{X, v) exists. 

Let denote the elapsed service time of the nth customer at time zeros and define = l^i 
if he is no longer in the system at time zero, then < Vn- Therefore E^^R{X^iy) < oo and 
E^^R{X,iy) <oc. 

As \{Rn{X + h,jy) - Rn{X,jy))/h\ < max^^^^^^<^<o K/A^ and |(i?^(A,z/ + h) - Rn{X,v))/h\ < 
max;^^ ^+^<n<o by Lebesgue Dominated Convergence Theorem, we have 

—E^R{X,u) = E^—R{X,u) and —E^R{X,u) = E^—R{X,u) 
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As the interarrival times have a continuous distribution, P(14,/z/ = |A^|/A) = for n < 0. 
Combining the change of hmit and the sample path analysis we have 

^E^R{\iy) = jE^E{X,iy) and ^E^R{\iy) = -^E^V{\iy) 

□ 

Table 3 shows the simulated results of an infinite server queue with base (i.e. A = 1) in- 
terarrival times distributed as Gamma(2, 2) and base (i.e. u = 1) service times distributed as 
Lognormal(— 0.25, 0.5). 



Table 3: Simulation result from exact sampling. 



(A, I.) 








f^E^R^{\v) 


(80, 1) 


7.0741 X 10-^ 


-1.1320 


6.1022 X 10-^ 


-2.8389 


(100, 1) 


5.6470 X 10-3 


-1.1316 


4.9379 X 10-3 


-2.9495 


(120, 1) 


4.7236 X 10-3 


-1.1337 


4.2337 X 10-3 


-3.0684 
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